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Summary 


)  A  horizontally  layered  system  of  fluids  and  solids  is 
bounded  above  and  below  by  half- spaces  which  may  be  solid  or 
fluid.  The  method  of  dynamic  stiffness  coupling  is  used  to 
obtain  a  matrix  relation  connecting  'spectral'  displacements 
and  external  stresses  at  the  interfaces.  The  pressure  due  to 
a  time-harmonic  acoustic  point  source  is  obtained,  via  this 
matrix  relation,  by  numerical  evaluation  of  a  Hankel  transform; 
'exactly'  using  a  Gaussian  quadrature  formula  and  approximately 
using  the  FFT  algorithm.  Numerical  results  demonstrate  that  the 
Fortran  programs  are  a  valuable  tool  for  predicting  short-range 


sound  propagation  at  low  frequencies. 
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INTRODUCTION 


An  important  and  challenging  research  area  in  the  field  of 
underwater  acoustics  is  the  problem  of  short-range  sound 
propagation  at  low  frequencies.  The  capability  to  predict  and 
understand  the  sound  characteristics  of  an  acoustic  noise  range 
would  facilitate  greatly  the  interpretation  of  the  measured  data 
which  are  invariably  contaminated  by  reflections  from  the  sea- 
surface  and  seabed;  additionally,  it  would  be  a  valuable  research 
tool  for  studying  the  variety  of  seismic  waves  that  may  occur  in 
a  layered  seabed. 

The  basic  mathematical  model  used  to  study  low-frequency 
propagation  is  the  Pekeris  model  Cl, 2, 31  which  consists  of  a 
time-harmonic  point  source  in  a  layer  of  fluid  which  is  bounded 
above  by  a  vacuum  and  bounded  below  by  a  half -space  of  fluid 
whose  density  and  sound  velocity  differ  from  those  of  the  fluid 
layer.  Clement  C43  has  given  some  numerical  results,  obtained 
from  a  Hankel  transform  by  numerical  integration,  of  the  Pekeris 
model,  and  has  also  included  for  comparison  purposes  numerical 
results  obtained  from  its  image  approximation,  which  is  a  simple 
closed-form  expression  obtained  by  summing  the  images  realized  by 
repeated  reflection  in  the  sea  surface  and  seabed.  However, 
because  the  seabed  is  likely  to  consist  of  a  layer  of  sediment 
overlaying  a  solid  substrate,  a  propagation  model  in  which  the 
seabed  is  multilayered  is  desirable. 

The  amount  of  previous  work  on  sound  propagation  in  multi¬ 
layered  media  is  too  great  to  reference  here.  Publications  by 
Mook  et  al.  C53  and  Schmidt  &  Jensen  C61  contain  a  comprehensive 
list  of  references.  In  most  multilayered  models  the  acoustic 
field  due  to  a  time  harmonic  point  source  excitation  is  evaluated 
by  numerical  integration  of  a  Fourier  or  Hankel  transform.  Mook 
et  al.,  in  an  interesting  paper  on  numerical  aspects,  demonstrate 
that  computation  is  much  more  reliable  and  quicker  if  the  poles 
associated  with  real  free-waves  in  fluid  layers  can  be  removed 
from  the  integrand  and  treated  explicitly:  however,  it  is  not 
clear  that  this  method  would  be  practical  when  elastic  layers 
which  include  dissipation  are  considered.  Schmidt  &  Jensen  have 
developed  an  efficient  method  for  solving  numerically  the  sound 
field  in  a  system  consisting  of  acoustic  fluid  and  elastic  solid 
layers:  they  use  the  fast  Fourier  transform  to  evaluate  the 
pressure  at  large  distances  from  the  source. 

The  approach  to  be  used  herein  is  the  method  of  dynamic 
stiffness  coupling  in  the  'spectral'  domain  which  has  been 
developed  at  ARE  (Teddington)  to  analyse  the  sound  radiation  from 
layered  media  excited  by  time-harmonic  point  forces  or  acoustic 
sources  C7,8,93.  The  method  generates  a  complex  frequency 
dependent  band-matrix,  which  relates  'spectral'  excitations  and 
displacements  at  the  layer  interfaces,  from  the  individual 
dynamic  stiffness  matrices  of  the  layer  elements.  Thus,  the 
method  has  much  in  common  with  the  well-known  finite  element 
method.  Moreover,  the  matrix  operations  required  to  give  the 
spectral  displacements  are  usually  better  conditioned  than  those 
of  the  transfer  matrix  method,  which  is  frequently  used  in 
problems  of  this  type. 


2.  PROBLEM  FORMULATION 


( a )  General 


The  mathematical  model  for  a  horizontally  stratified  ocean 
and  seabed  is  shown  in  Figure  1.  It  consists  of  an  arbitrary 
number  of  homogeneous  layers  whose  pressures  and  particle 
displacements  are  solutions  of  the  exact  linearized  equations  of 
acoustics,  elasto-dynamics  or  visco-dynamics.  It  is  excited  by  a 
time-harmonic  point  source  of  sound  which  is  located  in  one  of 
the  acoustic  fluids.  The  time-harmonic  factor  exp(-iwt),  where  uj 
is  the  radian  frequency,  is  omitted  from  all  equations. 


( b )  Hankel  Transforms 


For  the  axi-symmetric  problem  considered  here,  it  is 
convenient  to  represent  the  unknowns  by  Hankel  transforms  and 
their  inverses,  viz.. 
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In  the  above  equations  u„  and  u  are  the  radial  and  vertical 

L  Z 

components  of  the  particle  displacement,  and  p  is  the  acoustic 
pressure . 


( c )  System  Matrix  Equation 

The  spectral  dynamic  stiffness  matrices  of  the  various 
components  of  the  system  relate  'spectral'  stresses  or  pressures 
to  'spectral'  displacements  at  the  component's  interf ace ( s ) .  The 
matrices 


CSp(a)3,  CS£,(a)  3  ,  CSv(a)3  (2.3) 

2*2  4*4  4*4 

for  the  acoustic,  elastic  and  viscous  layers,  respectively,  are 
obtained  in  Sections  3-5.  The  matrices  for  the  upper  (U)  and 
lower  (L)  half-spaces, 

CSp( a) 3 ,  tsjj(a)3,  CS^(a)3  (2.4a) 

1*1  2*2  2*2 
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CSM(a) 3  EU(a) 3  =  CEM(a)3  (2.5) 

M*M  M*1  M*1 

where  M=2N+2  with  N  being  the  number  of  layers,  is  a  standard 
procedure  in  finite  element  theory  which  reflects  continuity  of 
displacement  and  equilibrium  of  stress  at  the  interfaces.  It  is 
discussed  in  Section  9.  The  matrix  CSM(a)3  is  the  system 
dynamic  stiffness  matrix  whose  elements  are  (in  general)  complex 
and  whose  bandwidth  is  7. 

In  equation  (2.5)  the  matrix  CEM(cx.)3  is  a  column  vector  of 
externally  applied  spectral  forces,  at  the  interfaces,  which  are 
due  to  the  time-harmonic  point  source  excitation.  Externally 
applied  stresses  are  defined  to  be  positive  when  acting  in  the 
positive  direction  of  the  coordinate  axes,  see  Figure  2.  This 
matrix  is  derived  in  Section  7  for  the  cases  in  which  the  source 
is  either  in  an  acoustic  layer  or  in  an  acoustic  half-space. 

The  spectral  displacements  at  the  N+l  interfaces 

T 

CU(a)  3  =  CUrl (a) ,Uzl (a) , . . . ,Ur(N+1 j (a) ,UZ(N+1 j (a) D  (2.6) 

are  found  by  matrix  inversion  of  equation  (2.5).  In  Section  8 
the  pressure  in  an  acoustic  layer  or  half-space  is  obtained  from 

knowledge  of  CU(oi)3. 


3.  ACOUSTIC  FLUID  LAYER 


( a )  General 

Figure  2A  shows  a  section  through  an  acoustic  fluid  layer 
with  lower  boundary  at  z=0  and  upper  boundary  at  z=h.  These 
surfaces  are  subject  to  prescribed  normal  'spectral'  pressures 

T 

Cp(a)3  =  Cp ( a ,h ) ,p ( a , 0 ) 3  (3.1) 

which  produce,  normal  to  the  surfaces,  the  'spectral' 
displacements 


T 

Cu  (oi)3  =  Cu  (a,h),u  (a, 0)3 

7  7  7 


A  matrix  CSp(a)3  is  required  which  relates  the  'spectral' 


(3.2) 


pressures  and  displacements  by  the  relation 


CSp(tx)  3Cuz(a)  3  =  Cp(<x)3 
2*2  2*1  2*1 


(3.3) 


< b )  Equations  of  Motion 

The  spectral  solution  of  the  linearized  wave  equation 


V2p(r,z)  +  k2p(r,z>  =  0 


(3.4) 


is  obtained  via  the  Hankel  transform,  equation  (2.1),  as 


p(a,z>  =  ( a ) exp( i\pZ )  +  A^ (a)exp( -iy^z )  (3.5) 

2  2  1/2 

where  "Yp  =  <  k  -a  )  with  lm(-yp)^0  in  order  to  satisfy  the 

radiation  condition;  k=u/c  is  the  acoustic  wavenumber;  u  is  the 
radian  frequency  and  c  is  the  sound  velocity  of  the  fluid  whose 
density  is  p. 

Equation  (3.5)  and  the  acoustic  momentum  equation 


evaluated  at 


the 


3p(a,z )  /3z  = 
layer  boundaries 


2“  , 

pw  u  ( a , z ) 

gives  the  matrix  equation 


(3.6) 


and 


"p  ( a ,  h ) ' 

_p  ( a ,  0  )  _ 

exp( iyph ) 
1 


exp( -iYph  > 
1 


"A  ^  ( a  >  ' 
.A 2  (a). 


(3.7) 


'u  ( a ,  h )  * 
z 

,  .  .  2 . 

rexp(  i-Tph) 

-exp(  -ivph)‘ 

[A  (a)' 

=  ( i\„/ pw  ) 

! 

u  (a,0 ) ' 
z 

r 

L  1 

-1 

[A  2  ( a ) . 

T 

from  which  CA^  ( a. )  , A2  (a )  3  may  be  eliminated  to  give  the  matrix 
equation  (3.3)  as 


2 

(pw  /-YpSin-Yph) 


'-cosYpb  1 

'u  (  a ,  h ) n 
z 

'p(a,h)‘ 

-1  cos-yph. 

u  ( a ,  0  ) 

L  z  J 

_p  (  a ,  0  ) . 

(3.9) 


u.a  j 
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4.  ELASTIC  LAYER 


( a )  General 

Figure  3A  shows  a  section  through  an  elastic  layer  with 
lower  boundary  at  z  =  0  and  upper  boundary  at  z=h.  These 
surfaces  are  subject  to  prescribed  tangential  and  normal 
'spectral'  stresses 

_____  T 

Ct(<x)3  =  Cx  (a,h),T  <a,h),T  ( oc, 0 ) ,t  (a,0 ) 3  (4.1) 

it  z  4  z  r  z  z  z 

which  produce  tangential  and  normal  'spectral'  displacements 

T 

Cu(a)3  =  Cu  ( a,h) ,u„ ( a ,h> ,u  ( a, 0 ) ,u  ( a, 0 ) 3  (4.2) 

r  z  r  z 


A  matrix  CSg(a.)3  is  required  which  relates  the  'spectral' 
stresses  and  displacements  by  the  relation 

C  S£ (  a, )  3 1  u  ( a )  3  =  Cx(a)3  (4.3) 

4*4  4*1  4*1 


( b )  Equations  of  Motion 

The  linear  elastic  equation  of  motion  C103 


u+m>v(v.u>  +•  \jv2u  =  pa2u/at2 


(4.4) 


(4.5) 


may  be  reduced  to  the  two  wave  equations 


V  F(r,z>  +  k£F(r,z)  =  0 


V2G(r,z)  +  k^,G(r,z)  =  0 


by  means  of  the  substitutions 


u  (r,z)  =  3F(r,z)/3r  +  a^Gt r ,z ) /3r3z 
u„(r,z)  =  3F(r,z)/3z  -  32G(r,z)/3r2  -  ( 1 / r ) 3G( r ,z ) / 3r  (4.6) 


ue=  0 


The  'spectral'  solutions  of  these  wave  equations  are 
obtained  via  the  zero  order  Hankel  transform,  equation  (2.1),  as 


F(a,z)  =  A^  (a)exp(  i-y^z  ) 


+ 


( a )  exp  (  -  i -y^z  ) 


9 


(4.7) 


G(ot,z)  =  A^  ( oc )  exp  (  i-y^z  )  +  A^  (  a )  exp( -iy^z  ) 

In  the  above  equations  x  and  m  are  the  Lame  constants  of 
he  elastic  material  whose  density  is  p ;  k^  is  the  wavenumber 

f  longitudinal  waves  in  the  elastic  solid;  k,p  is  the  wavenumber 

f  shear  waves  in  the  elastic  solid; 


=  to/c^  =  w/C (x+2u) /p3‘ 


kT  =  w/cT  =  w/(p/p) 


(4.8) 


■here  and  c^,  are  respectively  the  velocities  of 

ongitudinal  and  shear  waves;  (  k^-a2  )  ^  /  ^  and  =  (  k2-a2  )  ^ 1  2  , 

’ith  lm("y^)^0  and  lm(-y^)>>.0  in  order  to  satisfy  the  radiation 
ondit ions . 


c  >  Spectral  Displacements 

The  'spectral'  displacements  are  obtained  by  substituting 
he  Hankel  transform  representations  of  F  and  G  into 
■quations  (4.6),  and  noting  the  Bessel  function  identities 


x2d2J  (x)/dx2  +  xdJ  (x)/dx  +  (x2-n2)J  (x)  =  0 
n  n  n 


(4.9) 


dJQ ( x ) /dx  =  - ( x ) 


(4.10) 


biev  are 


u  (a,z)  =  -aF  -  a9G/3z 


u  ( a ,  z  )  =  SF ' 3z  +  a  G 
z 


(4.11) 


Ug ( a , z )  =  0 


and  G  may  be  eliminated  from  equations  (4.11)  by  the  use  of 
quations  (4.7)  to  give 


u  <a,z)  =  -txA1  (  a)exp(  i>r  z  )  -  aA_  (  a)  exp  ( -iyr  z  ) 

i  i.  Li  Z  L* 

-  ia-YrpA^  ( a. )  exp (  iyvpZ  )  +  iot>TA4  (  a )  exp( -i>Tz 


(4.12) 
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uz(ot,z)  =  ( a )  exp(  iy^z  )  -  1-yl”2  f  ex^ (  ~  ^L2 


+  a^A?  (a)  exp(  iYTz  )  +  a  A4  (  oc) exp( -iY^z  ) 


( d )  Spectral  Stresses 

It  is  straightforward  to  show  that  the  stresses 

T  =  y(3u  / 8 z  +  3u  /3r) 
rz  r  z 

t  =  X(  3u  /3r  +  u  /r  +  3u  /3z)  +  2u3u  / 3z 
zz  r  r  z  z 


(4.13) 


have  the  following  'spectral'  representations  in  terms  of  the 

unknowns  A.  to  A. . 

1  4 


x  =  -2iapYr  (oc)exp(  iYr  z  )  +  2iayYr  A_  (  a)  exp( -i Yr  z  ) 

if  Z  L  1  Lj  Li  Z  i_i 

2  2  °  2 
+  pafYrp-a.  )  A^  ( a.)  exp  (  iYrpZ  )  +  yafY^-a.  )  A4  ( ot )  exp( -iY,j.z  ) 

(4.14) 

—  2  2  ?  2 

=  -  (  Xk£+2yY'L )  A^  (a)  exp  (  iY^z  )  -  <  xk£  +  2yYL  )  A2  ( tx )  exp(  ~i^Lz  ) 

2  7 
+  2i<x  PYTA3(<x)exp(  iYTz)  -  2iot“yYTA4  (  a )  exp  ( -iYTz  ) 


( e )  Matrix  Equation 

Equations  (4.12)  and  (4.14),  evaluated  at  z=h  and  z=0, 
give  the  matrix  equations 


CR^.(a)3  CA(cx)3  =  Cu(a)3 
4*4  4*1  4*1 


[PjpfaJD  C  A  ( a )  3  =  C  x ( a) 3 
4*4  4*1  4*1 


(4.15) 


from  which  CA'a)3  may  be  eliminated  to  give  the  required 
relation  between  surface  'spectral'  stresses  and  displacements. 


viz.. 


CP_(a)3  CRc,(a)3  ^Cu(a)3  =  tx(ot)3 
E  E 


(4.16) 


The  elements  of  the  matrices  C P-^. f  ot )  3  and  CPg(a)3  are  listed 
the  Appendix. 


VISCOUS  FLUID  LAYER 


1  General 


Figure  3A  shows  a  section  through  a  viscous  fluid  layer  with 
«?er  boundary  at  z  =  0  and  upper  boundary  at  z=h.  The  surfaces 
?  subject  to  prescribed  tangential  and  normal  'spectral' 

resses  Cx(a)3  which  produce  'spectral'  surface  displacements 

(a)].  A  matrix  CS^(a)D  is  required  which  relates  the  surface 

pectral'  stresses  and  displacements  by  the  equation 

CSy(a)]  Cu(al]  =  [x(a)]  (5.1) 

4*4  4*1  4*1 


)  Equations  of  Motion 

The  linearized  Navier-Stokes  equations  of  a  viscous  fluid  CIO] 
-Vp  +  yV2u  +  (y/3)V(V.u)  =  p3u/3t 

3p/3t  +  pV.u  =  0  (5.2) 

3  p / 3t  =  C23p/3t 

y  be  reduced  to  the  equations 

V2F(r,z)  +  Cw2/ ( c2-4iuy/3p) DF< r ,z >  =  0 

V2G(r,z)  +  ( iwp/y )G( r ,z )  =  0  (5.3) 

p  =  - ( ipc2/w) ( V . u ) 


means  of  the  substitutions 


u  (r,z>  =  3F(r,z)/3r  +  3  G(r,z)/3r3z 
r 

u  (r,z)  =  3F(r,z)/3z  -  32G(r,z)/3r':'  -  (  1  /  r  )  3G(  r  ,  z  )  /  3r 


(5.4) 


In  the  above  equations  y  is  the  dynamic  coefficient  of 
3cosity  of  the  fluid  whose  density  is  p,  and  in  which  the 

und  speed  is  c.  The  overdot  (')  denotes  differentiation 
th  respect  to  time. 

Equations  (5.3)  and  (5.4)  together  with  the  stress-velocity 
la t ions 


x  =y(3u/3z+3u/3r) 
rz  r  z 


(5.5) 


.  1 I  (  *  I  If 


C  the  sound  levels  exceed  the  datum.  The  plots  are  consistent 
h  two  free-waves  in  the  water  whose  amplitudes  decrease  with 
tance  from  the  free-surface  and  the  sediment  interface, 
pectively . 

The  situation  is  more  complex  at  5,  7  and  9Hz  ,  shown  in 
■ures  7-9.  Here  the  increasing  number  of  wavenumber  branches 
:es  physical  interpretation  difficult  without  the  availability 
eigenvectors.  The  most  interesting  feature  in  these  plots  is 
•  short-wavelength  interference  pattern  obtained  when  both  the 
irce  and  hydrophone  are  near  to  the  bottom:  this  is  due  to  the 
eraction  between  the  short-wavelength  free-wave  which  is 
■dominantly  a  shear  wave  in  the  sediment  (labelled  1  in  Figure 
,  and  the  other  waves  whose  wavelengths  are  much  larger. 


Ef feet  of  Shear 


Figures  10  and  11  show  the  sound  level  versus  distance  plots 
ained  by  Gaussian  integration  when  the  sediment  layer  and 
ie r  half-space  are  both  treated  as  acoustic  fluids.  The 
>quencies  chosen  for  these  plots  were  5  and  9Hz ,  respectively, 
romparison  of  Figures  7  and  10,  which  are  plots  for  5Hz,  shows 
it  the  effect  of  shear  is  strongest  when  both  the  source  and 
!  receiver  are  near  to  the  bottom:  as  expected,  the  short- 
'elength  interference  effect  disappears  when  shear  is  omitted, 
romparison  of  Figures  9  and  11,  which  are  plots  for  9Hz ,  shows 
ior  differences  at  all  three  source/ receiver  positions,  the 
>t  striking  being  in  the  case  in  which  both  source  and  receiver 
>  near  to  the  bottom. 


Sound  Levels  by  FFT  Integration 

The  fast  Fourier  transform  algorithm  has  also  been  used  to 
iluate  the  Hankel  transform.  Plots  corresponding  to  those  in 
fures  5-9  have  been  obtained  but  they  are  not  reproduced  here 
:ause  they  are  almost  indistinguishable  from  those  obtained  by 
issian  quadrature  at  3,  5,  7  and  9Hz;  and  at  1Hz  they  differ  by 
to  3dB  -  which  result  is  not  surprising  because  the  FFT 
suits  are  strictly  valid  only  when  kr>l. 

The  Gaussian  method  was  found  to  be  more  efficient  than  the 
?  method  because  whereas  the  former  requires  a  high  sampling 
isity  of  the  integrand  only  near  to  f ree-wavenumber s ,  the 
;ter  requires  a  high  (equally  spaced)  sampling  density  through- 
;  the  range  of  integration.  One  of  the  disadvantages  of  the 
issian  method  is  the  large  number  of  Bessel  function  evaluat- 
is  required,  but  this  could  be  somewhat  alleviated  by  inter- 
.ation  or  computations  of  restricted  accuracy.  A  subsequent 
slication  will  give  a  detailed  comparison  of  numerical  results 
:ained  by  the  two  methods,  and  will  also  discuss  the  parameters 
:essary  for  numerical  integration. 
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propagating  free-waves  to  exist,  the  values  of  and  y^,  in 

the  underlying  elastic  half -space  must  he  purely  imaginary.  Thus, 
the  real  wavenumbers  of  free-waves  are  constrained  to  lie  in  the 
region 


< rock )) k^< rock )  (11.1) 

and  the  wave  amplitudes  in  the  lower  half-space  must  decay 
exponentially  with  distance  frun  the  interface. 

Each  branch  of  the  dispersion  plot  represents  free-waves 
whien  are  various  combinations  of  longitudinal  waves  in  the  water 
layer  and  shear  and  longitudinal  waves  in  the  sediment  layer  and 
rock  half-space.  The  relative  magnitudes  of  these  waves  may  be 
determined  exactly  by  calculating  the  eigenvectors  of  the  matrix 
CSM(a)3  for  each  pair  of  values  of  a  and  u,  but  some  general 
speculations  are  possible  without  doing  so. 

As  the  frequency  is  increased  the  character  of  each  free- 
wave  changes.  The  branch  labelled  1  starts  as  a  coupled  wave 
whose  motion  is  predominantly  longitudinal  in  the  sediment  layer 
and  rock  half-space,  it  then  changes  rapidly  to  a  predominantly 
shear  wave  in  the  sediment.  The  branches  labelled  2  and  3  start 
as  predominantly  longitudinal  waves  in  the  sediment,  turn  into 
acoustic  waves  in  the  fluid,  and  then  later  change  into 
predominantly  shear  waves  in  the  sediment.  The  branch  labelled  4 
starts  as  a  longitudinal  wave  in  the  sediment  and  turns  into  a 
shear  wave  in  the  sediment.  The  branch  labelled  5  starts  as 
predominantly  a  shear  wave  in  the  sediment  and  turns  into  a 
longitudinal  wave  in  the  sediment. 

Figure  4B  shows  the  real  branches  of  the  radial  wavenumber 
versus  frequency  plot  for  the  system  in  which  shear  is  absent  in 
the  sediment  layer  and  rock  half -space,  i.e.  they  are  both 
treated  as  fluids.  The  wavenumbers  are  constrained  to  lie  in  the 
region 


a>k( rock)  (11.2) 

and  they  are  not  too  dissimilar  from  the  wavenumbers  that  would 
have  been  obtained  in  the  absence  of  the  fluid  sediment  layer.  A 
comparison  of  Figures  4A  and  4B  shows  that  shear  plays  a  signifi¬ 
cant  role  in  determining  the  wavebearing  properties  of  the 
system . 


( c )  Sound  Levels  by  Gaussian  Integration 

Figures  5-9  show  the  sound  level  versus  distance  plots 
obtained  by  Gaussian  integration  for  the  three  (A,B,C)  source/ 
receiver  positions.  At  1Hz,  shown  in  Figure  5,  the  sound  level 
decreases  rapidly  with  distance  but  tends  to  level  off  when  the 
hydrophone  is  located  on  the  bottom.  This  is  consistent  with  the 
presence  of  a  coupled  free-wave  in  the  sediment  layer  and  rock 
half -space,  which  has  only  a  weak  component  in  the  water  layer. 

At  3Hz,  shown  in  Figure  6,  the  sound  levels  in  Plot  A  are 
much  less  than  those  of  the  spherical  spreading  datum:  in  Plots  B 
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.  NUMERICAL  RESULTS 


)  General 


The  material  and  geometric  constants  in  SI  units,  used  in 
e  computation  of  sound  propagation,  in  a  layer  of  water  on  a 
abed  consisting  of  a  layer  of  sediment  on  top  of  a  hard  sub- 
ratum,  were  chosen  as  follows: 


>per  Half -Space:  vacuum 

.yer  of  Water:  p  =  1000 

c=1500 

h=225 

iyer  of  Sediment:  p=1400 

c=cl=1750  ct= 

=  250 

X=4. llxlC9 

y=8 . 75x10  7 

>wer  Half -Space:  p=2500 

c=cL=5000  ct= 

=  2500 

Rock)  A=3.13x1010 

p= 1 . 56x10 1 0 

Distance  off  Seabed 

Plot  A 

Plot  B 

Plot 

Source 

125 

125 

10 

Hydrophone 

100 

1 

1 

h=20 


ius  in  the  'A'  plot  of  Figures  5-11  both  source  and  hydrophone 
re  near  to  midwater;  in  the  'B'  plot  the  source  is  near  midwater 
id  the  hydrophone  is  near  to  the  bottom;  and  in  the  'C  plot 
Dth  the  source  and  the  hydrophone  are  near  to  the  bottom.  The 
nplitude  of  the  point  source  was  chosen  as  Pq=1 ;  hence  the 

Dund  levels  in  dB  ref.  1  micropascal  at  the  stated  range  can  be 
sferred  to  a  free-field  source  strength  of  120dB  at  lm.  The 
iund  levels  are  shown  as  the  source  moves  along  a  horizontal 
rack  ranging  from  x=0  (point  of  closest  approach  to  hydro- 
ione)  to  x=1000m,  the  horizontal  separation,  y,  between  source 
id  hydrophone  at  the  point  of  closest  approach  being  selected  as 
10m.  Each  of  the  plots  contains  the  spherically  spreading  case, 

=20.01og10< l/RJ+120,  where  R=CxZ+y Z+< z-zQ > Z 3 1 7 2  for  comparative 

jrposes.  Damping  was  included  by  setting  the  hysteretic  loss- 

actors  to  0.01,  a  value  which  does  not  introduce  unduly  high 

ttenuations  at  the  maximum  range  plotted.  The  infinite  limit  of 

itegration  in  the  Hankel  transform  was  truncated  to  a.  =0.40 

max 

lich  is  sufficiently  high  to  include  all  the  f ree-wavenumbers . 

iwever,  due  to  the  limited  maximum  number  size  ( 1 0 3 s )  on  the 
3P-11  computer  it  was  necessary  to  subdivide  the  water  layer 
ito  three  to  prevent  numerical  overflow  caused  by  the  term 
up ( -iyph) . 


b )  Wavenumber  Plots 


Figure  4A  shows  the  real  branches  of  the 
ersus  frequency  plots  of  the  dissipationless 
ree-waves  cut-on  at  0.0,  2.2,  3.5,  6.5  and  9 


radial  wavenumber 
system,  in  which 
0Hz.  For  radially 


Here,  the  real  values  of  a  such  that  Re ( detCSMf a > 3 ) 
vanishes  are  found  by  a  simple  algorithm  which  searches  for  sign 
changes.  If  Im( detCSM( a) D >  is  also  zero  for  the  computed  value 
of  <x  then  the  wavenumber  of  a  free-wave  has  been  found. 


( c )  Gaussian  Quadrature 

The  infinite  limit  in  the  integral,  equation  (10.1),  is 

truncated  to  a  finite  value,  a  ,  which  is  determined  either 

max 

from  a  wavenumber  versus  frequency  plot  or  by  trial  and  error. 

The  integral  is  then  evaluated  by  an  adaptive  Gaussian  quadrature 
(order  2)  scheme,  which  splits  the  range  of  integration  into  a 
user  selected  number  of  equal  intervals  each  of  which  is 
repeatedly  halved  until  a  relative  convergence  test  is  satisfied 
by  successive  approximations  to  the  integral  in  the  interval. 
Large  savings  in  computation  time  are  possible  by  parallel 
computation  at  an  array  of  r-values,  due  to  the  simplicity  of  the 

JQ((xr)  term  in  relation  to  p(oc,z). 


<d>  Fast  Fourier  Transform 


The  integral  in  equation  (10.2)  is  proportional  to 

00 

I(r>  =  |  f (a)exp( iar )da 
-00 


(10.4) 


1/2- 

where  f(a)=a  p(a.z).  The  discrete  version  of  this  equation 
has  been  given  by  Cooley  et  al.  C12D  in  a  form  which  is  suitable 
for  evaluation  by  a  fast  Fourier  transform  algorithm,  viz.. 


N-l 


I  (  r  ) 


=  Sa  y  p^.exp<  2iri  ik/N) 


k=0 


Pk 


00 


f ( kSa+NqSa) 
-00 


(10.5) 


M 

where  N  =  2  with  integer  M;  r^=j5r  for  j  =  0  to  N-l;  and 

S<x=2tt /NSr .  For  pressure  evaluation  the  computation  time  is 
halved  when  the  relation 


f ( -a )  =  if ( a)  (10.6) 

is  used  when  equation  (10.5)  is  evaluated.  The  procedure  for 

choosing  the  parameters  q  ,  N,  Sa  (and  hence  Sr)  is  to  be 

max 

discussed  elsewhere. 
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The  reduced  matrices  which  occur  in  the  case  of  fluid  layers 
and  half-spaces,  and,  if  a=0,  also  in  the  case  of  elastic  and 
viscous  layers  and  half-spaces,  require  no  special  attention 
before  assembly  into  the  system  stiffness  matrix.  However,  a 
singular  system  matrix  which  contains  all  zeros  in  certain  rows 
may  be  generated,  in  which  case  the  relevant  diagonal  terms  must 
be  set  to  unity. 


10.  COMPUTER  IMPLEMENTATION 


(a)  Programs  for  Acoustic  Pressure 

Fortran  programs  have  been  written  to  evaluate  numerically 
the  sound  pressure  level  in  an  acoustic  layer /half -space  due  to 
the  presence  of  a  time  harmonic  source  located  in  an  acoustic 
layer/half -space.  One  program  evaluates  the  pressure  by  Gaussian 
quadrature  of  its  Hankel  transform  representation,  viz., 

00 

p(r,z)  =  <  1  /  2  tt  )  |  p(a,z)  JQ((xr  )ada  (10.1) 

0 

and  another  program  evaluates  the  pressure  from  its  asymptotic 
expansion  for  large  r,  viz., 

00 

p(r,z)  ~  ( l/8ir3r  )  1/2exp< -iir/4)  j  p(  a,  z  )  exp(  iar  )  ot1  /  2dtx  (10.2) 

-OO 

by  the  use  of  the  fast  Fourier  transform  algorithm. 

Because  the  transform  representations  are  arbitrary  with 
respect  to  a  solution  of  the  corresponding  homogeneous  problem, 
it  is  necessary  to  introduce  damping  into  the  problem  in  order  to 
lift  any  singularities  off  the  real  axis,  thus  permitting 
numerical  integration.  This  is  done  by  setting 

X=X  (  1-iri  )  ,  p=p(  1-in  )  ,  c  =  c(l-ir)  )  (10.3) 

A  C 

where  the  n-values  are  the  hysteretic  loss-factors. 


( b )  Free-Wave  Propagation 

In  order  to  aid  the  physical  interpretation  of  numerical 
results  a  Fortran  program  has  been  written  to  calculate 
wavenumber  versus  frequency  dispersion  curves.  In  the  absence  of 
external  sources  and  dissipation,  the  system  of  homogeneous 
equations  (2.5)  has  a  non-trivial  solution  if  and  only  if 
detCSM(a)]  vanishes.  For  a  given  value  of  u>  there  will,  in 
general,  be  both  real  and  complex  values  of  a  for  which  the 
determinant  vanishes.  Real  values  of  a  are  the  wavenumbers  at 

1/2 

which  free-waves  propagate,  cylindrically  decaying  as  1/r 
Complex  values  of  a  may  represent  evanescent  waves  whose 
amplitudes  decrease  exponentially  with  distance. 


CSE(a)U 

4*4 


( pur/YpSinYph) 


0 

0 

0 

0 


0 

cosvph 

0 

1 


0 

0 

0 

0 


0 

-1 

0 

•cosYph 


(9.3) 


( b )  Upper  and  Lower  Half-space  Elements 

It  is  also  convenient,  before  assembly,  to  write  the 
acoustic  half-space  matrices  in  the  same  form  as  the  elastic  and 
viscous  half-space  matrices.  Thus 


and 


CSp(a)  3Cu(  a)  3  =  Cp(a)3 
1*1  1*1  1*1 


CSp(  ot)  3Cu(  a)  3  =  Cp(<x)3 
1*1  1*1  1*1 


(9.4) 


(9.5) 


are  expanded  into  the  form 


and 


CSp(  a)  3Cu(  a)  3  =  Cx(oc)3 
2*2  2*1  2*1 


(9.6) 


in  which 


CS^(a)3Cu(a)3  =  Ct(<x>3 
E 


2*2  2*1 


2*1 


(9.7) 


and 


CSp( a ) 3  =  <ipw2/YF>  [  o  J  ] 
CS^(a)3  =  ( -ipu2 / Yp)  [  °  J  ] 


(9.8) 


(9.9) 


( c )  Computer  Implementation 


The  assembly  of  the  component  matrices  to  form  the  system 
dynamic  stiffness  relation  is  a  standard  finite  element  procedure 
which  reflects  the  continuity  of  displacement  and  the  equilibrium 
of  stress  at  the  component  interfaces.  Computer  implementation 
is  straightforward  for  devotees  of  the  finite  element  method. 

The  sign  convention  that  externally  applied  stresses  are  defined 
to  be  positive  when  acting  in  the  positive  direction  of  the  co¬ 
ordinate  axes  makes  it  necessary  to  change  the  sign  of  the 

elements  of  the  last  two  rows  of  each  layer  matrix,  CSE<a>3  and 

CS..(a)3,  and  to  change  the  sign  of  all  the  elements  of  the  upper 
v  — n  — u 

half-space  matrix,  CSp(a)3  or  CS^(a)3,  before  assembling  the 

component  matrices.  An  example  of  an  assembled  system  matrix  is 
shown  elsewhere  C73. 
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—  2  — 

p(<x,z)  =  -i(pw  /-Yp)uz(a,0)exp(  iypZ) 

+  (  2iripQ/-Yp)  Cexp(  ivpl  2-z0  |  )  +  exp( iyp( z+zQ ) ) D  (8.2) 

For  the  case  in  which  there  is  no  acoustic  source  in  the  upper 
half-space,  it  is  necessary  only  to  set  P0=0  in  equation  (8.2). 


( c )  In  Lower  Half-space 

The  spectral  displacement  uz(a,0)  of  the  boundary  of  the 

lower  half-space  is  obtained  from  the  dynamic  stiffness  relation, 
equation  (2.5).  If  the  lower  half -space  contains  a  source,  the 
spectral  pressure  is  obtained  from  equations  (7.11)  and  (7.12) 
evaluated  at  z=0,  the  boundary  of  the  half -space.  Thus  after 
some  algebra. 


—  2  — 

p(a,z)  =  i(pw  /-yp)uz(«.,0)exp( -ivpZ ) 

+  (  2iripQ/Vp)  Cexp(  i>-pl  z-zQ  |  )  +  exp( -i\p(  z+zQ  )  )  3  (8.3) 

For  the  case  in  which  there  is  no  acoustic  source  in  the  lower 
half-space,  it  is  necessary  only  to  set  p0  =  0  in  equation  (8.3). 


9.  ASSEMBLY  OF  ELEMENTS 


( a )  Laver 

Figure  2  shows  the  standard  sign  convention  for  stresses  and 
pressures.  Both  displacements  an<*  externally  applied  stresses 
are  defined  to  be  positive  when  acting  in  the  positive  direction 
of  the  coordinate  axes.  Note  that  the  directions  of  positive 
stress  and  positive  pressure  are  opposite. 

It  is  convenient,  before  assembly,  to  write  the  acoustic 
layer  matrices  in  the  same  form  as  the  elastic  and  viscous  layer 
matrices.  Thus 

CSp( a) 3Cu( a) 3  =  Cp(a)3  (9.1) 

2*2  2*1  2*1 

is  expanded  into  the  form  of  equation  (4.3) 

CS£(a)  DCu(  a.)  D  =  Ct(<%)3  (9.2) 

4*4  4*1  4*1 


in  which 
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Evaluation  of  equations  (7.11)  and  (7.12)  at  z  =  0,  the 
boundary  of  the  lower  half -space,  gives  the  relation 


CS„(a)3Cu  (a)  3  =  Cp(a)3  +  CE^(a)3  (7.13) 

ir  Z  r 

in  which 

CEp3  =  -4iripQexp( -ivFz0> /-yf  (7.14) 

Equations  (7.13)  and  (7.14),  together  with  the  sign 
convention  which  requires  external  stresses  to  be  positive  when 
acting  in  the  positive  direction  of  the  coordinate  axes,  show 
chat  the  point  source  in  the  lower  half-space  is  equivalent  to  an 

external  vertical  stress  of  -Ep  at  the  boundary  of  the  lower 

half -space . 


8.  PRESSURE  IN  ACOUSTIC  FLUIDS 


(a)  In  Fluid  Laver 

The  spectral  displacements,  u  (a,h)  and  u  (<x,0),  of  the 

£•  4- 

layer  boundaries  are  obtained  from  the  dynamic  stiffness  matrix 
relation,  equation  (2.5).  If  the  layer  contains  a  source,  the 
spectral  pressure  is  obtained  from  equations  (7.3)  and  (7.4) 
evaluated  at  the  layer  boundaries,  z  =  h  and  z=0.  Thus  after 
some  considerable  algebra 

—  2  —  — 

p(cx,z)  =  (  pw  /-ypSinyph)  Cuz  ( a,0 )  cos\p(  h-z  )  -  u„  ( a,h)  cos-ypZ3 

-  (  27rp0/-YpSinYph)  CcosYp(h- j  z-z0  |  )  +  cosypt h-z-zQ ) 3  (8.1) 

For  the  case  in  which  there  is  no  acoustic  source  in  the  layer, 
it  is  necessary  only  to  set  Pq  =  0  in  equation  (8.1). 


(b)  In  Upper  Half-space 

The  spectral  displacement  u  (<x,0)  of  the  boundary  of  the 

upper  half-space  is  obtained  from  the  dynamic  stiffness  relation, 
equation  (2.5).  If  the  half-space  contains  a  source,  the 
spectral  pressure  is  obtained  from  equations  (7.7)  and  (7.8) 
evaluated  at  z=0,  the  boundary  of  the  half -space.  Thus  after 
some  algebra 


Ep2  =  4irp0cos(YFCh-z03) /vFsin(Vph) 

Equations  (7.5)  and  (7.6),  together  with  the  sign  convention 
which  requires  external  stresses  to  be  positive  when  acting  in 
the  positive  direction  of  the  coordinate  axes,  show  that  the 
point  source  in  a  fluid  layer  is  equivalent  to  external  normal 
stresses  of  -Epi  and  Ep2  at  the  upper  and  lower  boundaries 

respectively. 


( c )  Source  in  Upper  Fluid  Half-space 

The  'spectral'  form  of  equation  (7.2)  in  the  upper  half- 
space  is  obtained  from  equations  (7.1),  (7.2)  and  (6.1)  as 


p(a,z)  =  2iripQexp(  ivFl  z-zQ  |  ) />F  +  A^ ( a) exp( iypZ )  (7.7) 

and  the  'spectral'  form  of  the  momentum  relation  is 
2— 

pu>  uz(cx,,z)  =  -2irp0sgn(  z-zQ  )  exp<  i\Fl  z-z0  |  )  +  iypA.^  ( <x)  exp(  i^pZ  ) 

(7.8) 

Evaluation  of  equations  (7.7)  and  (7.8)  at  z=0,  the 
boundary  of  the  upper  half-3pace,  gives  the  relation 


CSp(a)DCuz(a)3  =  Cp(a)3  +  CE^(ot)3  (7.9) 

in  which 

CEp3  =  -4iripQexp(  iYFzQ ) /yf  (7.10) 

Equations  (7.9)  and  (7.10),  together  with  the  sign 
convention  which  requires  external  stresses  to  be  positive  when 
acting  in  the  positive  direction  of  the  coordinate  axes,  show 
that  the  point  source  in  the  upper  half-space  is  equivalent  to  an 

external  vertical  stress  of  Ep  at  the  boundary  of  the  upper 

half -space . 


(d )  Source  in  Lower  Fluid  Half-space 

The  'spectral'  form  of  equation  (7.2)  in  the  lower  half- 
space  is  obtained  from  equations  (7.1),  (7.2)  and  (6.4)  as 

p(ot,z)  =  2iripQexp(  iyF  j  z-z0  |  ) /yF  +  A2  (a)exp( -i-ypZ)  (7.11) 

and  the  'spectral'  form  of  the  momentum  relation  is 
2— 

pw  uz<a,z)  =  -2irpQsgn(  z-zQ  )exp(  iYp!  z-zQ  |  )  -  iYpA2  ( a) exp( -ivpZ  ) 

(7.12) 
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7.  EXCITATION  STRESS  OF  POINT  SOURCE 


7.  EXCITATION  STRESS  OF  POINT  SOURCE 

• 

(a)  Green's  Function 

The  pressure  of  a  point  source,  located  in  acoustic  fluid 
at  ( r  ,  z )  =  ( 0  ,zQ) ,  see  Figure  2,  is  best  developed  via  the 

• 

' 

transform  of  its  free-space  Green's  function  CUD. 

*  •  *  \  * 

00 

• 

—  « 

pQexp(  ikRQ) /RQ  =  ipQ  J  ( 1/ yf>  Jq ( otr  )  exp(  i*yF |  z-zQ  |  )<xdtx 

0 

(7.1) 

which  is  augmented  by  a  scattering  term  when  boundaries  are 
present.  Thus  the  total  pressure  is 

• 

— 

p( r ,z )  =  pQexp(iJcR0) /Rq  +  ps(r,z) 

(7.2) 

-  . 

where  pg(r,z)  satisfies  the  homogeneous  reduced  wave  equation 
(3.1).  In  the  above  p^  is  the  amplitude  of  the  point  source 

W 

,  • 

- 

2  2  1/2 

and  RQ=Cr  +<z_2o)  is  the  distance  from  the  source  to 

the 

observation  point. 

(b)  Source  in  Fluid  Layer 

l 

The  'spectral'  form  of  equation  (7.2)  is  obtained  from 
equations  (7.1),  (7.2)  and  (3.5)  as 

•v 

p(a,z)  =  2irip0exp(  iypl  z-zQ  |  ) /yF 

-  ‘v  ’• 

+  A-^  ( a)  exp(  i-YpZ  )  +  A2  ( a) exp( -i\pZ  ) 

(7.3) 

"i 

and  the  'spectral'  form  of  the  momentum  relation,  equation 
is 

(3.3) , 

• 

2— 

pw  u2<<x,z)  =  -27rp0sgn(  z-zQ  )exp(  ivp|  z-z0  |  ) 

.. 

*.  .'I 

+  ifpCA^  (oc)exp(  i>pZ  )  -  A2  ( a)  exp( -i>pZ  )  3 

(7.4) 

t 

-• 

Evaluation  of  equations  (7.3)  and  (7.4)  at  the  layer 
boundaries,  z=h  and  z=0,  gives  the  matrix  relation 

•»! 

CSp( a) 3  Cuz(a)3  =  Cp(a)3  +  CEp(a)3 

(7.5) 

- ! 

2*2  2*1  2*1  2*1 

V*\  - 

in  which  the  elements  of  CEp(a)3  are 

*.< 

Epi  =  4irp0cos(7pZ0)/>'pSin(-Yph) 

(7.6) 

_  • 

‘  1 

17 


=  2iauv^A?  ( a )  exp(  -iv^:  )  +  pa'\^,-a  )  (  a  )  exp  (  -i>^,z  ) 

_  2  2  2 

t__  =  -  ( \k^+2u>^ )  A^  ( a  >  exp( -iv^z  )  -  2ia  u-y^lalexpl  -iv^z: 


(6.14) 


Equations  (6.13)  and  (6.14),  evaluated  at  2=0,  the  boundary  of 
the  lower  halt-space  give  the  matrix  equations 

CR^( a) 3CA( a) 3  =  Cu(a)3 


(6.15) 


CP^(  a)  3CA<  a)  3  =  Ct(cx>3 


from  which  CA(a)3  may  be  eliminated  to  give  the  required 
relation  between  surface  'spectral'  stresses  and  displacements 


CP^(a)3CR^(a)3  XCu(a)3  =  Cx(a)3 
E  E 


(6.16) 


The  elements  of  the  matrices  CP^(a)3  and  CR^(a)3  are  listed 
the  Aooendix. 


(  f )  Viscous  Half-spaces 

The  method  of  deriving  the  matrices  CSy(<x>3  and  CS^(cx.)3 

is  the  same  as  that  described  in  Section  5.  Equations  (6.7)  and 
(6.12),  the  spectral  solutions  of  the  reduced  wave  equations  (5.3 
in  the  upper  and  lower  half-spaces  respectively,  are  substituted 
into  equations  (5.4)  and  (5.5).  Subject  to  the  identities  (5.8) 
the  spectral  dynamic  stiffness  matrices  of  the  viscous  half¬ 
spaces  are 


CS^( a ) 3  = 


-iwCP^(cx)3CR^(<x)3_1 


CS,L7(<x>3  =  -iwCP?t<  a)  a)  3  1 


(6.17) 


(6.18) 


in  the  upper  half -space. 


<  a, z ) 
uz (a, z ) 


-aA^ (a) exp( iy^z )  -  ia 

2 

i^L^i  (oc>  exp(  i^L2  *  +  a  A^  ( a) exp(  iy,pZ  ) 


yTA3<a)exp( iy,pZ ) 


(6.8) 


and 


2  2 


xr_  =  -2iayy^A1  ( a) exp(  iy^z  )  +  pcx(YT-oi.  ) A?  ( a) exp(  iy,p2  ) 

—  2  2  2 

tzz  =  -(\kL+2yyL)A1((X)exp(iyLz)  +  2ia  yy^A^  (ot)exp(  iyTz ) 


(6.9) 


Equations  (6.8)  and  (6.9)  evaluated  at  z  =  0,  the  boundary  of  the 
upper  half -space,  give  the  matrix  equations 


CR^(  a.)  JCA(a)  3  =  Cu(a)D 


2*2  2*1 


,U, 


2*1 


(6.10) 


CPr,(a)DCA(oc)D  =  Cx(tx)D 
E 


2*2  2*1 


2*1 


from  which  CA(a)3  may  be  eliminated  to  give  the  required 
relation  between  surface  'spectral'  stresses  and  displacements, 
viz .  , 


CPg(a)3ERg(a>3  1Cu(cx)D  =  Ct(oc)3 


(6.11) 


2*2  2*2 


2*1 


2*1 


The  elements  of  the  matrices  CP^(a)3  and  CRg(<x)3  are  listed 
the  Appendix. 


( e )  Elastic  Lower  Half-space 

The  spectral  solutions  of  the  reduced  wave  equations  (4.5) 
representing  outgoing  waves  in  the  lower  half-space  are 


F(a,z)  =  A2 (a) exp( -iy^z ) 
G(a,z)  =  A^ (a.) exp( -iy^z  ) 


(6.12) 


which,  when  substituted  into  equations  (4.11)  and  (4.13),  give 
the  following  expressions  for  the  spectral  displacements  and 
stresses  in  the  upper  half-space. 


ur<<x,z)  =  -<xA2  (<x)exp( -iyLz>  +  ia,yTA4  ( a)  exp< -iyTz  ) 

—  2 

uz(a,z)  =  -iyLA2  (a)  exp! -iy^z  )  +  a  A^  ( a)  exp( -iy^,z  ) 


(6.13) 
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the  spectral  stresses  and  displacements  for  both  upper  and  lower 
elastic  and  viscous  half-spaces,  respectively. 


(b)  Acoustic  Fluid  Upper  Half-space 

The  spectral  solution  of  the  linearized  wave  equation  (3.4) 
representing  outgoing  waves  in  the  upper  half-space  is 

p(a,z)  =  (ot)  exp(  i\pZ  )  (6.1) 

Equation  (6.1),  together  with  the  acoustic  momentum  equation 
(3.6),  evaluated  at  the  boundary  of  the  half-space,  z=0,  allows 
the  spectral  pressure  to  be  related  to  the  spectral  displacement 
of  the  boundary,  viz.. 


where 


ESp( a) 3Cu( a) 3  =  Cp(a,0)3 
1*1  1*1  1*1 


Sp( a)  = 


■ipu)  / Yp 


(6.2) 


(6.3) 


( c )  Acoustic  Fluid  Lower  Half-space 

The  spectral  solution  of  the  linearized  wave  equation  (3.4) 
representing  outgoing  waves  in  the  lower  half-space  is 

p(a,z)  =  A2(a)exp( -iYpZ)  (6.4) 


Equation  (6.4),  together  with  the  acoustic  momentum  equation 
(3.6),  evaluated  at  the  boundary  of  the  lower  half -space,  z=0, 
allows  the  spectral  pressure  to  be  related  to  the  spectral 
displacement  of  the  boundary,  viz., 

CSp<  a)  3Cu(  cx)  3  =  Cp(a,0)3  (6.5) 

1*1  1*1  lvi 

where 

Sp(a)  =  ipw2/YF  (6.6) 


(d)  Elastic  Upper  Half-space 


The  spectral  solutions  of  the  reduced  wave  equations 
representing  outgoing  waves  in  the  upper  half-space  are 

F(a,z)  =  A^ ( a) exp( iy^z ) 

G(a,z)  =  A3 (a) exp( iyTz ) 


(4.5) 


(6.7) 


which,  when  substituted  into  equations  (4.11)  and  (4.12)  give  the 
following  expressions  for  the  spectral  displacements  and  stresses 
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-t  =  X  ( 2  3u  /  3z  +  /r  +  30  /  3z)  +  2y30  /3z 
zz  r  r  z  z 

enable  the  analysis  to  proceed  as  in  Section  4  to  give  the  matrix 
relation 

CPy(ot)3  CA(a)3  =  Ct(<x>3 
4*4  4*1  4*1 


(5.6) 


CRy(a)D  CA(a)3  =  Cu(a)3 
4*4  4*1  4*1 

from  which  CA(a)3  may  be  eliminated  to  give  the  required 
relation  between  'spectral'  stresses  and  displacements,  viz.. 


-iwCPy( a) 3  CRy(a) 3  1  Cu(cx)3  =  Cx(a)3  (5.7) 
4*4  4*4  4*1  4*1 

The  elements  of  the  matrices  CPy(a)3  and  CRv(cx>3  are  the  same 
as  those  of  CPg,(a.)3  and  CRp(a)3  for  the  elastic  layer,  except 
that 


=  w2/ (c2-4iwy/3p)  -  a2 

2  2 
Yt  =  (iwp/p)  -  a. 

2 

x  =  ipc  /u  -  2y/3 


(5.8) 


M  =  U 


6.  THE  UPPER  AND  LOWER  HALF- SPACES 


( a )  General 

j  Figure  2B/2C  shows  a  section  through  an  acoustic  half-space 

whose  boundary,  z=0,  is  subject  to  prescribed  normal  'spectral' 

'  pressure  which  produces  a  'spectral'  displacement  normal  to  the 

!  surface.  Matrices  CSp(o03  and  CSp(ot)3  are  required  which 

:  1*1  1*1 

relate  the  spectral  pressures  and  displacements  in  the  upper  and 
’  lower  half-spaces.  Figure  3B/3C  shows  a  section  through  an 

elastic  or  viscous  half-space  whose  boundary,  z=0,  is  subject  to 
*■  prescribed  tangential  and  normal  'spectral'  stresses  which 

produce  tangential  and  normal  'spectral'  displacements.  Matrices 

I  _  CSg(a)3,  CSp(a)D,  CSy(a)3  and  CSy(a)3  are  required  which  relate 

2*2  2*2  2*2  2*2 


( 


4 


4 


4 


« 


« 


< 


_  4 


4 
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12.  CONCLUDING  REMARKS 


The  mathematics  of  acoustic  propagation  in  layered  media  has 
been  presented  in  a  form  which  is  not  too  difficult  to  implement 
on  a  minicomputer  by  programmers  familiar  with  the  finite  element 
method.  In  addition  to  illustrating  the  importance  of  including 
shear  in  the  seabed  parameters,  the  limited  number  of  numerical 
results  presented  here  have  demonstrated  the  capability  of  the 
Fortran  program  as  a  potentially  useful  tool  for  predicting  low 
frequency  sound  propagation  characteristics  at  short  ranges. 

Some  additional  numerical  work  along  the  following  lines  is 
proposed  as  follow-up  projects: 

(i)  comparison  of  results  obtained  from  various 
integration  algorithms  including  those  used 
herein.  In  particular,  numerical  problems 
concerned  with  explicit  treatment  of  poles  C5D 
should  be  investigated  when  both  shear  and 
dissipation  are  present. 

(ii)  numerical  studies  of  acoustic  radiation  using 
geometric  and  material  constants  typical  of 
acoustic  measuring  ranges,  and  the  applicability 
of  wavenumber  versus  frequency  plots  to  the 
physical  interpretation  of  the  results. 

^iii)  computation  of  the  particle  displacements  in  both 
the  sea  and  seabed,  which  can  be  done  by  the 
existing  program  with  trivial  modification.  The 
displacements  in  the  sea  enable  acoustic  intensity 
vectors  to  be  calculated  C3D,  and  the  displace¬ 
ments  in  the  sediment  would  facilitate  the  study 
of  low  frequency  seismic  interface  waves. 


Of  special  value  would  be  the  availability  of  suitable  experi¬ 
mental  data  obtained  from  noise  ranges  or  scale  model  tests, 
which  would  help  to  determine  the  range  of  applicability  of  the 
theoretical  work  contained  herein. 


E  A  Skelton  (HSO) 


Manuscript  completed  January  1985 
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FIG. 5  SOUND  LEVEL  VERSUS  DISTANCE  (X) 
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FIG. 6  SOUND  LEVEL  VERSUS  DISTANCE  (X) 
FREQUENCY=3HZ .  WITH  SHEAR. 


FIG. 7  SOUND  LEVEL  VERSUS  DISTANCE  (X) 
FREQUENCY* 5HZ .  WITH  SHEAR . 
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FIG. 10  SOUND  LEVEL  VERSUS  DISTANCE  (X) 
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APPENDIX 


The 

Elastic  Layer  Matrices 

CPE(a) 3 

and  CRtr(a)3 

E 

• 

P 

Ell 

=  -2iyayEexp( iy^h) 

p 

E12 

=  2iyayEexp( -iy^h) 

PE1 3 

2  2 

=  ya(  Yrp-a.  )exp(iyTh) 

p 

rE14 

2  2 

=  ya(yT-a  )exp(-iyTh) 

• 

PE21 

=  - ( Xk^+2yyE ) exp( iy^h) 

PE2  2 

=  - ( Xk^+2uyE >exp( -iy^h) 

PE23 

2 

=  2iya  yTexp(iyTh) 

p 

E24 

2 

=  -2iytx  yTexp( -iyTh) 

• 

Rows  3  and  4  are  obtained  by  setting  h=0  in  rows  1  and  2. 


RE11 

=  -aexpfiy^h) 

RE1 2 

=  -txexp(  -iy^h) 

PE13 

=  -ia.yTexp(  iyTh) 

RE14 

=  iocyTexp( -iyTh) 

RE21 

=  iyLexp(iyLh) 

RE22 

=  -iyLexp( -iy^h) 

RE2  3 

2 

=  oc  exp(  iyrph ) 

RE24 

2 

=  a,  exp( -iyTh) 

Rows  3  and  4  are  obtained  by  setting  h=0  in  rows  1  and  2. 

For  the  case  when  a=0  the  matrix  CRE(cx)3  is  singular  and 
matrix  inversion  is  not  possible.  However,  the  non-zero  elements 

of  CSg(0)3  can  be  shown  to  be  -* 


^E22  =  (  *+2y )  k^cosk^h/sink^h 
=  -  (  X+2y  )  k^/ sinkjJn 
Sg42  =  ( x+2y ) k^/ sink^h 
^E44  =  - ( X+^M ) k^cosk^h/ sink^h 
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The  Elastic  Half- space  Matrices 


CPg(a)D 


and 


CFg(ot)D 


,U  T . 

Eli  =  -2i^>l 

PU 

El  2 

=  ya(\2-a2) 

E21  =  -(XkL+2^L5 

PU 

E22 

_  .  2 
=  2iy<x  vT 

:E11  '  - 

RU 

El  2 

=  -ia\T 

-  iv 

E21  i 'L 

RU 

*E22 

-  oc2 

For  the  case  when  ot=0  the  matrix  CRg(<x)D  is  singular  and 
matrix  inversion  is  not  possible.  However,  the  non-zero  element 

of  CSg(0)3  can  be  shown  to  be 

®E22  ■  1<x+2»)kL 


The  Elastic  Half-space  Matrices 


tp|(a)3 


and  CR^(a)D 


Ell  =  2iM<xvL 

pk 

rEl  2 

=  ya(y^,-a2) 

E21  *  -UkL+2^L’ 

pk 

E22 

=  -2iyot2YT 

Ell  = 

rl 

E12 

—  i  ocyij* 

'E21  =  ~lyL 

rl 

*E22 

=  a2 

For  the  case  when  a=0  the  matrix  CR^(a)D  is  singular  and 
matrix  inversion  is  not  possible.  However,  the  non-zero  element 

of  ESg(0)U  can  be  shown  to  be 


-i( A+2y)kr 
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